/*
 *  elasticlog.cpp
 *  EPPI-0.0
 *
 *  Created by Sergio Conde Martin on 04/05/11.
 *  Copyright 2007 __MyCompanyName__. All rights reserved.
 *
 */

#include "elasticlog3d.h"

ElasticLog3D::ElasticLog3D()
:	HyperelasticStretches(),
	E(0.0),
	nu(0.3),
	lambda(0.0),
	mu(0.0)
{}

// ***

ElasticLog3D::ElasticLog3D(ifstream &in, const string& na, const string & t)
:	HyperelasticStretches(na,t),
	E(0.0),
	nu(0.0),
    mu(0.0),
    lambda(0.0),
    K(0.0),
    b(2.0),
    typeVol("")
{
	string keyword;
	while (in >> keyword && keyword.compare("ENDMATERIAL") )
	{
		if(!keyword.compare("VOL"))
		{
			in >> typeVol; 
			if(!typeVol.compare("mixed")) in >> b;
		}
		if(!keyword.compare("E"))	in >> E;
		else if(!keyword.compare("NU"))	in >> nu;
        else if(!keyword.compare("MU"))	in >> mu;
        else if(!keyword.compare("LAMBDA"))	in >> lambda;
        else if(!keyword.compare("K"))	in >> K;
		else if(!keyword.compare("DENSITY"))	in >> density;
	}
    if (E != 0.0 and nu != 0.0)
    {
        lambda = nu * E / ((1.0 + nu) * (1.0 - 2.0 * nu));
        mu = E / ( 2.0 * (1.0 + nu) );
        //K = E / (3.0 * (1.0 - 2.0 * nu));
    }
    if (K == 0.0)
    {
        K = lambda + 2.0/3.0 * mu;
    }
	
	hypertype = "coupled";
	
	if(!typeVol.compare("cuad") || !typeVol.compare("log") || !typeVol.compare("mixed"))
	{
		hypertype = "decoupled";
	}
	if (!hypertype.compare("decoupled"))
	{
		Svol.setZero(6);
		Siso.setZero(6);
		Ctvol.setZero(6,6);
		Ctiso.setZero(6,6);
	}
}


// ***

void ElasticLog3D::strainEnergy()
{
	if(J > 0.0)
	{
		if(!hypertype.compare("coupled"))
		{
			W = mu * (log(lambda1) * log(lambda1) + log(lambda2) * log(lambda2) + log(lambda3) * log(lambda3)) + lambda / 2.0 * log(J) * log(J);
		}
		else
		{
			Wiso = mu * (log(lambda1m) * log(lambda1m) + log(lambda2m) * log(lambda2m) + log(lambda3m) * log(lambda3m));
			
			if(!typeVol.compare("cuad"))	Wvol = K / 2.0 * (J - 1.0) * (J - 1.0);
			if(!typeVol.compare("log"))	Wvol = K / 2.0 * log(J) * log(J);
			if(!typeVol.compare("mixed"))	Wvol = K / (b * b) * (1.0 / pow(J,b) - 1.0 + b *log(J));
			
			W = Wvol + Wiso;
		}
	}
	else
	{
		cout << "*ERROR: Negative volume, material " << name << " has collapsed \n";
		W = -1.0;
	}
}

void ElasticLog3D::DstrainEnergy()
{
	if(J > 0.0) 
	{
		if(!hypertype.compare("coupled"))
		{
			W1 = 2.0 * mu * log(lambda1) / lambda1 + lambda * log(J) / lambda1;
			W2 = 2.0 * mu * log(lambda2) / lambda2 + lambda * log(J) / lambda2;
			W3 = 2.0 * mu * log(lambda3) / lambda3 + lambda * log(J) / lambda3;
		}
		else
		{
			Wiso1 = 2.0 * mu * log(lambda1m) / lambda1m;
			Wiso2 = 2.0 * mu * log(lambda2m) / lambda2m;
			Wiso3 = 2.0 * mu * log(lambda3m) / lambda3m;
			
			if(!typeVol.compare("cuad"))	dWvol = K * (J - 1.0);
			if(!typeVol.compare("log"))	dWvol = K * log(J) / J;
			if(!typeVol.compare("mixed"))	dWvol = K / (b * J) * (1.0 - 1.0 / pow(J,b));
		}

	}
	else
	{
		W1 = - 1.0;
		W2 = - 1.0;
		W3 = - 1.0;
	}
}

// ***

void ElasticLog3D::DDstrainEnergy()
{
	if(J > 0.0)
	{
		if(!hypertype.compare("coupled"))
		{
			W11 = ( 2.0 * mu * (1.0 - log(lambda1)) + lambda * (1.0 - log(J)) ) / ( lambda1 * lambda1 );
			W12 = lambda / (lambda1 * lambda2);
			W13 = lambda / (lambda1 * lambda3);
			W22 = ( 2.0 * mu * (1.0 - log(lambda2)) + lambda * (1.0 - log(J)) ) / ( lambda2 * lambda2 );
			W23 = lambda / (lambda2 * lambda3);
			W33 = ( 2.0 * mu * (1.0 - log(lambda3)) + lambda * (1.0 - log(J)) ) / ( lambda3 * lambda3 );
		}
		else
		{
			Wiso11 = 2.0 * mu * (1.0 - log(lambda3m));
			Wiso12 = 0.0;
			Wiso13 = 0.0;
			Wiso22 = 2.0 * mu * (1.0 - log(lambda2m));
			Wiso23 = 0.0;
			Wiso33 = 2.0 * mu * (1.0 - log(lambda3m));
			
			if(!typeVol.compare("cuad"))	ddWvol = K; 
			if(!typeVol.compare("log"))	ddWvol = K * (1.0 - log(J)) / (J * J);
			if(!typeVol.compare("mixed"))	ddWvol = K / (b * J * J) * ( (b + 1.0) / pow(J,b) - 1.0);
		}
	
	}
	else
	{
		W11 = - 1.0;
		W12 = - 1.0;
		W13 = - 1.0;
		W22 = - 1.0;
		W23 = - 1.0;
		W33 = - 1.0;
	}
}

// ***

void ElasticLog3D::DDDstrainEnergy()
{}

// ***

double ElasticLog3D::getE()
{	return E;}

// ***

double ElasticLog3D::getNu()
{	return nu;}

// ***

double ElasticLog3D::getLambda()
{	return lambda;}

// ***

double ElasticLog3D::getMu()
{	return mu;}

// ***

double ElasticLog3D::getDensity()
{	return density;}

// ***

double ElasticLog3D::getBulkmodulus()
{	return K;}

// ***
